function [LaserWL, ramanpeak, LaserWLgauss] = paracetamolcalibration(datasetname, SIF)

% [LaserWL] = paracetamolcalibration(datsetname, SIF)
%
% This function uses the spectrum of Paracetamol (acetamidophenol) to 
% return the exact wavelength of the exitation Laser to be used for 
% subsequent calibration of actual Raman spectra.
%
% INPUT:    datasetname     =   Name of the C.H.-spectra file
%           SIF             =   if set 1 input will be treated as ANDOR
%                                  *.sif file and parset to SIFREAD first
%
% OUTPUT:   LaserWL         =   Wavelength of exitation Laser 
%           ramanpeak       =   Raman peak chosen for calibration by user

if nargin<1,
    help paracetamolcalibration
    return
end

%% -----==================== CONFIG ====================-----
% ----- Set the reference peaks -----
ramanpeak0 = 797.2;         % put here in cm-1
ramanpeak1 = 857.9;         % put here in cm-1
ramanpeak2 = 1168.5;         % put here in cm-1
ramanpeak3 = 1236.8;         % put here in cm-1
ramanpeak4 = 1278.5;         % put here in cm-1
ramanpeak5 = 1323.9;         % put here in cm-1
ramanpeak6 = 1371.5;         % put here in cm-1
ramanpeak7 = 1561.6;         % put here in cm-1
ramanpeak8 = 1648.4;         % put here in cm-1
ramanpeak = ramanpeak5;

% Reference spectrum
refx = [650 ...
        680.8 700.8 710.8 720.8 740.8 ...
        767.2 787.2 797.2 807.2 820 ...
        820 827.9 834.5 847.9 857.9 867.9 887.9 ...
        1139 1159 1169 1179 1199 ...
        1206.8 1226.8 1236.8 1246.8 1260 1264.8 ...
        1268.5 1278.5 1288.5 1300 ...
        1300.1 1313.5 1323.5 1333.5 ...
        1341.5 1361.5 1371.5 1381.5 1401.5 ...
        1495.1 1510.1 1515.1 1520.1 1535.1 ...
        1541.1 1555.1 1561.1 1568.1 1586.1 ...
        1600 1620 1630 1635 ...
        1639.4 1648.4 1658.4 1678.4 ...
        ];
    
refy = [0 ...
        0 3 10 3 0 ...
        0 5 22 5 0 ...
        0 3 5 3 44 5 0 ...
        0 5 38 5 0 ...
        0 5 40 3 8 ...
        3 22 5 1 3 ...
        1 5 58 5 0 ...
        5 20 5 0 ...
        0 1 3 1 0 ...
        0 5 20 5 2 ...
        3 50 5 3 ...
        4 45 5 0 ...
        ];

% sets point to fit on Gaussian    
pointsongaussianEdit = 50;


%% -----==================== CODE ====================-----
ramanpeak = ramanpeak/10000000;

%% Check for identity of input format
if nargin==2
    if SIF == 1
        rc = loadsif_nonGUI(datasetname, 1);
        if rc ~= 0
            errordlg('''SIF''could not be read!!!')
            return
        end
        disp(['Sif read with the RC ',num2str(rc)])
    else
        errordlg('Variable ''SIF'' was not set correctly!!')
        return
    end
end


%% load *.mat data
%disp('start load mat')
try load(datasetname)
catch ME
  errordlg('ERROR: file not found in "cyclohexcalibration"!')
  return
end
load(datasetname)

% ----- pick peak -----
scrsz = get(0,'ScreenSize');
chooser = figure('OuterPosition',...
                [scrsz(3)*0.05 (scrsz(4)*0.05) ...
                    scrsz(3)*0.9 scrsz(4)*0.9]);
    subplot(2,1,1);
    plot(datax,datay); 
%    h = zoom;   
%    zoom
%     q = msgbox({'Zoom in on the most prominent peak, please...' '' '' ...
%                 'Press OK when the peak maximum is easy to click.'},...
%                 'Zoom in on Peak');
%     posq = get(q,'Position'); posq(1,2) = 100;
%         set(q, 'Position', posq);
%     uiwait(q)
%    setAllowAxesZoom(h,refgraph,false);
    q = msgbox('Click on the center of the peak, please...',...
                'Choose Peak', 'modal');
        posq = get(q,'Position'); posq(1,2) = 100;
            set(q, 'Position', posq);
    uiwait(q)

    [x y] = ginput(1);
    text(x, y,'selected peak','Color','r','EdgeColor','r', ...
        'FontWeight', 'bold','HorizontalAlignment','center');
    zoom('out')
    
 %% MARTINS Gaussfit

    gc=x;                     % Generate short dataset for Gaussian fit
    vFitwidth=floor(pointsongaussianEdit/2);
    gci=find(datax>gc, 1);
    gk=datax(gci-vFitwidth:gci+vFitwidth);
    gd=datay(gci-vFitwidth:gci+vFitwidth);
    ftype = fittype('gauss1');
    foptions = fitoptions('gauss1');
%    foptions.StartPoint = [2e4 gc 0.2];
    gfit = fit(gk,gd,ftype,foptions)
    gf=gfit.a1*exp(-((gk-gfit.b1)/gfit.c1).^2);
   %   gf=(gfit.a1*exp(-((gk-gfit.b1)/gfit.c1).^2))+(gfit.a2*exp(-((gk-gfit.b2)/gfit.c2).^2));
    %gf=10000*exp(-((gk-gc)/.1).^2);

   % Display the fit
   line(gk,gf,'Color','r','Marker','+','Linewidth',2,'Linestyle','none');
       
   

   
    
%% Plot Reference in second subplot
    subplot(2,1,2);
        plot(refx, refy, 'r')
         text(refx(2),refy(2),['\leftarrow',ramanpeak],...
      'HorizontalAlignment','left')
%   setAxesZoomMotion(h,refgraph,'horizontal');

%% anotate Referencespectrum:
% Create textarrow
annotation(chooser,'textarrow',[0.255729271791962 0.255729271791962],...
    [0.283704966665861 0.244598821414464],'TextEdgeColor','none',...
    'String',{'797.2 cm^-^1'});

annotation(chooser,'textarrow',[0.286831812255541 0.294654498044329],...
    [0.398268398268398 0.366522366522367],'TextEdgeColor','none',...
    'String',{'858 cm^-^1'});

annotation(chooser,'textarrow',[0.479791395045632 0.493363809440681],...
    [0.373737373737374 0.333178150217257],'TextEdgeColor','none',...
    'String',{'1168.5 cm^-^1'});

annotation(chooser,'textarrow',[0.529986962190352 0.539113428943937],...
    [0.38961038961039 0.346320346320346],'TextEdgeColor','none',...
    'String',{'1237 cm^-^1'});

annotation(chooser,'textarrow',[0.566492829204694 0.561929595827901],...
    [0.29004329004329 0.243867243867244],'TextEdgeColor','none',...
    'String',{'1278.5 cm^-^1'});

annotation(chooser,'textarrow',[0.61864406779661 0.601043024771838],...
    [0.402597402597403 0.431457431457431],'TextEdgeColor','none',...
    'String',{'1323.9 cm^-^1'});

annotation(chooser,'textarrow',[0.643415906127771 0.631473272490221],...
    [0.262626262626263 0.235411255411255],'TextEdgeColor','none',...
    'String',{'1371.5 cm^-^1'});

annotation(chooser,'textarrow',[0.827249022164276 0.807692307692308],...
    [0.408369408369408 0.375180375180375],'TextEdgeColor','none',...
    'String',{'1648.4 cm^-^1'});

annotation(chooser,'textarrow',[0.750977835723598 0.750977835723598],...
    [0.294372294372294 0.232323232323232],'TextEdgeColor','none',...
    'String',{'1561.6 cm^-^1'});

%% Ask for identity of selected peak

    choicelist = {[num2str(ramanpeak0),' cm^-^1'], ...
                    [num2str(ramanpeak1),' cm^-^1'], ...
                    [num2str(ramanpeak2),' cm^-^1'], ...
                    [num2str(ramanpeak3),' cm^-^1'], ...
                    [num2str(ramanpeak4),' cm^-^1'], ...
                    [num2str(ramanpeak5),' cm^-^1'], ...
                    [num2str(ramanpeak6),' cm^-^1'], ...
                    [num2str(ramanpeak7),' cm^-^1'], ...
                    [num2str(ramanpeak8),' cm^-^1']};
    [choice, ~] = listdlg('PromptString','Which peak did you choose?',...
                'Name','Choose Peak', ...
                'SelectionMode','single', ...
                'ListSize', [150 100], ...
                'InitialValue', 6, ...
                'ListString', choicelist);

%% Handle response
    switch choice
        case 1
            ramanpeak = ramanpeak0;
        case 2
            ramanpeak = ramanpeak1;
        case 3
            ramanpeak = ramanpeak2;
        case 4  
            ramanpeak = ramanpeak3;
        case 5
            ramanpeak = ramanpeak4;
        case 6
            ramanpeak = ramanpeak5;
        case 7
            ramanpeak = ramanpeak6;
        case 8
            ramanpeak = ramanpeak7;
        case 9
            ramanpeak = ramanpeak8;
        otherwise 
            errordlg('Callibration failed !!')
    end            

%% close Figures and compute WL
close(chooser)

% ----- Compute LaserWL -----    
LaserWL = (1/((ramanpeak/10000000)+(1/x)));
LaserWLgauss = (1/((ramanpeak/10000000)+(1/gfit.b1)));

%% output Message Dialog with results
results = {['Choosen Raman peak was: ',num2str(ramanpeak),'cm^-^1'],...
            '', ...
            ['Resulting LASER wavelength: ', num2str(LaserWL),' nm'],...
            '',...
            ['Resulting LASER wavelength (Gauss): ',num2str(LaserWLgauss),' nm']};
 q = warndlg(results,...
             'Results', 'modal');
 uiwait(q)